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Abstract. We present a simple physically motivated picture for the mildly non-linear regime 
of structure formation, which captures the effects of the bulk flows. We apply this picture 
to develop a method to significantly reduce the sample variance in cosmological N-body 
simulations at the scales relevant to the Baryon Acoustic Oscillations (BAO). The results 
presented in this paper will allow for a speed-up of an order of magnitude (or more) in 
the scanning of the cosmological parameter space using N-body simulations for studies which 
require a good handle of the mildly non-linear regime, such as those targeting the BAO. Using 
this physical picture we develop a simple formula, which allows for the rapid calculation of the 
mildly non-linear matter power spectrum to percent level accuracy, and for robust estimation 
of the BAO scale. 



1 Introduction 



The mildly non-linear regime of structure formation has grown in importance recently since 
the Baryon Acoustic Oscillations (BAO) in the matter power spectrum lie precisely in this 
regime at low redshift {z < 3), when dark energy comes to dominate. Ever since their 
detection by the SDSS [1], the BAO have become one of the central targets of cosmological 
investigations. 

The BAO arise (e.g. [2]) from the sound waves present in the tightly-coupled pho- 
ton/baryon fluid in the early universe. After recombination, the radiation pressure can no 
longer support the baryons, which until then have been driven away from overdensities. Thus, 
the baryons quickly lose momentum. The resulting comoving sound horizon corresponds to 
the scale of the BAO, which is ~ 150 Mpc. The BAO in the matter two-point function appear 
as a peak at that scale with a width of around 20 Mpc. That peak translates to a pattern of 
decaying oscillations in the matter power spectrum with a wavelength of about 0.06/i/Mpc 
(see Figure 1). 

This acoustic signature acts as a standard ruler which may allow us to probe the behavior 
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Figure 1: Shown schematically are various matter power spectra at 2; = for ACDM. The power 
spectra are divided by a smooth BBKS [3] power spectrum with shape parameter F = 0.15 in order 
to highlight the wiggles due to the BAO. The non-linear power spectrum (i.e. the "exact" power 
spectrum obtained from N-body simulations) is given by Pnl', the linear power spectrum by Pl; the 
part due to the "memory of the initial conditions" by R^Pl] the power spectrum calculated in the 
Zel'dovich approximation by P^; the power due to the projection of the non-linear density field on the 
Zel'dovich density field by B?Pz- The green (long dashed) arrow represents the power generated from 
mode-coupling due to bulk flows, free-streaming and structure formation; while the blue (dot-dashed) 
arrow represents the power generated by structure formation alone. The current Hubble expansion 
rate in units of lOOkm/s/Mpc is given by h. 
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of dark energy. Any evolution of the equation of state of dark energy may be detectable by 
a shift in the BAO scale. Current and future experiments such as BOSS^ and WFIRST^ 
are expected to reach 1% errors in the matter power spectrum at the scales relevant for the 
BAO; and sub-percent precision in the acoustic scale. 

The acoustic scale is well in the linear regime. Thus, one can expect that the position of 
the acoustic peak in the 2-point function will be only slightly affected by small-scale physics. 
However, reaching sub-percent-level accuracy in the acoustic scale at low redshift requires 
a good handle of the mildly non-linear regime, which (in Fourier space) is probed by the 
successive peaks and troughs in the power spectrum (Figure 1). One well-recognized effect 
from the mildly non-linear regime is the suppression of the BAO in the power spectrum due to 
matter flows, free streaming and non-linear evolution [2]. These effects introduce a ~ lOMpc 
broadening of the acoustic peak, which in turn degrades the accuracy in the determination 
of the acoustic scale (in real space). Understanding this degradation can help us construct 
better peak reconstruction methods, in analogy with [4, 5]. 

The effects of mode-coupling also introduce a percent-level bias in the location of the 
acoustic peak at low redshift [6], exactly where dark energy has an impact. Thus, in order to 
utilize the coming observations we need to extend our standard cosmological model at least 
to this mildly non-linear regime (scales of ~ lOMpc at z = 0). 

Going beyond the linear regime is a non-trivial task. Structure formation is a problem 
without analog, made especially difficult by the free streaming of Cold Dark Matter (CDM) 
and the long range of the gravitational interaction. Indeed, for arbitrary initial configurations 
solving the non-linear structure of CDM is a problem much more complicated than turbu- 
lence, as the latter can be treated entirely in configuration space as a fluid, using the ffist 
several velocity moments of the one-particle distribution function (due to the collisions and 
short range of the fluid interactions). While a consistent understanding of CDM structure 
formation in principle requires one to work with the full one-particle distribution function in 
phase-space. 

The reason why structure formation is at all tractable (beyond simulations) at early 
times is because of the extremely low entropy of its initial conditions (almost homogeneous 
and isotropic, with almost Gaussian perturbations of the gravitational potential). Moreover, 
even at late times many interesting phenomena (such as the BAO) are at scales which are 
either linear or only mildly non-linear. Since particles in Eulerian space have not moved a 
lot (~ lOMpc on average), at large scales this causes the hierarchy for the velocity moments 
for the CDM to collapse and one is left with the usual continuity and Euler equations (e.g. 
[7]). 

That result can in principle be nullified in the presence of strong backreaction effects. 
However, recently the effects of the small non-linear scales on the large linear scales have been 
carefully analyzed analytically in [7], and it was shown that virialized structures completely 
decouple from large scale modes, while non-virialized small scales introduce only an effective 
pressure and viscosity to the large-scale equations of motion. Thus, even in the presence 
of non-linearities, the linear large-scale modes can still be represented as an almost perfect 
effective fluid. 

At the scales relevant for the BAO, Eulerian Standard Perturbation Theory (SPT) (e.g. 
[8]) starts breaking down, and thus alternatives must be devised. One way of addressing the 
problem is through numerical simulations. This has been possible thanks to the exponentially 

^http:/ /www. sdss3.org/surveys/boss.php 
^http:/ /wfirst. gsfc.nasa.gov/ 
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increasing available computational power and the development of successively more efficient 
algorithms to solve the N-body problem. Thus, conducting cosmological experiments through 
simulations has proliferated over the last two decades. 

However, computational time is still a very important constraint on simulations. Sam- 
pling variance requires one to perform numerous simulations of different realizations of the 
same cosmology; or alternatively simulate large volumes. This prevents an efficient sampling 
of the cosmological parameter space, and thus, it is still difficult to get a complete handle 
on the uncertainties that would arise from measurements at the mildly non-linear scales rel- 
evant to the BAO (e.g. [9]). Thus, devising analytical or semi-analytical schemes to study 
the mildly non-linear regime is of great importance. Moreover, a good physical picture of 
what happens in this transitional regime is still needed. That may allow us to construct more 
robust statistics to study this regime, which are especially needed because what we observe 
are in fact tracers of the underlying density field, and not the field itself. 

1.1 The Zel'dovich approximation 

A parallel route for studying the mildly non-linear regime is using analytical techniques (e.g. 
[10]). One of the first successful models beyond the fluid approach of Eulerian SPT was the 
Zel'dovich approximation [11], which is easily tractable analytically. In this approximation, 
the CDM particles are allowed to move on fixed straight trajectories with particle velocities 
which are drawn from an irrotational stochastic vector field. Particles are allowed to intersect 
each other's trajectories. After the particle streams start crossing, a rich structure in phase- 
space develops. 

We will argue that the Zel'dovich approximation (ZA) is a crucial stepping stone in 
understanding the mildly non-linear regime. It is the first order in Lagrangian Perturbation 
Theory (LPT) (e.g. [12, 13]) where the expansion parameter is usually taken to be the 
overdensity, which makes the solutions tractable. In Eulerian space, however, we will see that 
the ZA is equivalent to certain infinite partial resummations of higher order contributions, 
which encode some of the effects of the large-scale modes on the mildly non-linear regime. 

We demonstrate the superiority of the ZA over linear theory^ in Figure 2, where we 
show the cross-correlation coefficient Pnl,z between the non-linear overdensity^ and the 
overdensity calculated in the ZA; as well as the cross-correlation coefficient pnl,l between 
the non- linear overdensity and the overdensity in linear theory. As one can see, pnl,l decays 
quickly in the mildly non- linear regime; while pnl,z remains close to 1 (see also [14]). Given 
such a high cross-correlation between the non-linear density field and the density field in 
the ZA, one may ask whether the ZA (and higher order LPT) can be used to correct for 
the uncertainties associated with cosmic variance at linear and even mildly non-linear scales 
(corresponding to the BAO) in simulations. This question will be addressed in Section 3. 

Let us now try to understand why linear theory (unlike the ZA) fails to produce a large 
cross-correlation with the exact result at the mildly non-linear scales, despite the fact that it 
gives better agreement than the ZA in the power spectrum at these scales (Figure 1). The 
reason for that will turn out to be the large-scale flows, which are well modeled by the ZA, 
but not by linear theory. 

^By "linear theory" we refer to the first order SPT. 

*The cross-correlation coefficient between two fields Sa and 5b is defined as j = {3a5ij)^ /{{\5a\'^){\Sb\^)). 
^By non-linear overdensity and power spectrum we refer to the "true" quantities calculated using N-body 
simulations for example. 
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Figure 2: The cross-correlation coefficient between the non-linear overdensity and the overdensity 
in linear theory and the ZA as a function of scale at z = 0. The cross-correlation coefficient, Pnl,l^ 
is very close to the non-linear propagator in RPT, despite their different definitions. Note that the 
Zel'dovich density field is well-correlated with the non-linear density field well into the mildly non- 
linear regime. 



The effect of the large-scale flows results in corrections in SPT, which are regulated 
by the parameter kay{< k), where av{< k) gives the root mean square (rms) particle dis- 
placement integrated up to a wavevector k (see also Section 2).^ Note that ay is given by 
the integral over the velocity power spectrum and is proportional to the velocity dispersion, 
which is why for brevity we will usually refer to as the velocity dispersion. The corrections 
from the largest-scale flows {k — )• 0) due to ay{< k) cancel at each order in SPT [15]. This 
must be the case because uniform motions translate coherently the density field, which has 
no observable effect when equal-time statistics are considered. However, this is no longer 
true when one considers statistics at different times. 

Linear theory predicts a density field which is simply an overall rescaling of the initial 
density field. Therefore, pnl,l also equals the cross-correlation between the initial and final 

®This quantity is defined as (fc(T„(< fc))^ = kikj{viVj){< k) — (47r/3)fc'^ PL{k')dk' , where Pl is the linear 
matter power spectrum. This implies that is in fact equal to (l/\/3)x(rms particle displacements). 
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Figure 3: We show the overdensity 6{< k) and rms particle displacements CTt,(< k) integrated up 
to fc, calculated in linear theory at z = 0. The quantity fccr„(< k) is due to the bulk flows and 
it regulates the destruction of the cross-correlation between the linear and non-linear density fields. 
Note that ka^{< k) is larger than d{< k) in the mildly non-linear regime. The peak of the velocity 
power spectrum is at fc* which is well below the non-linear scale fc^vL when the linear density power 
per logarithmic fc-interval reaches 1. Therefore Uy is roughly constant in the mildly non- linear regime. 
In the plot we marked with vertical lines the positions of k^ and k^L f^or reference. See the text for 
further discussion. 

density fields. This is a statistics at different times and so it is affected by the bulk flows. 
The velocity power spectrum peaks at rather large scales (at k^. 0.07/i/Mpc; see Figure 3), 
therefore it is precisely large-scale motions that dominate. So, in the mildly non-linear 
regime the density field at a later time can be crudely thought of as a translated version of 
the initial density field, where the translation is a random variable. Such a uniform random 
shift destroys the cross-correlation between the initial and final density fields. This explains 
why pnl,l decays so quickly in the mildly non- linear regime. 

In Figure 3 we show the magnitude of kay[< k) and the overdensity, 5{< k), integrated 
up to k (the latter is defined through eq. (5.2)). Since the velocity power spectrum peaks 
at large scales, the velocity dispersion (t^(< k) reaches a constant well below the non-linear 
scale, kNL, (see Figure 3). Thus, even though kay{< k) grows large in the mildly non-linear 
regime, its effects should be accounted for by the ZA, which exactly captures the effects of 
the convective derivative in the Euler equation. In the mildly non- linear regime kav{< k) is 
larger than 5{< k). Therefore, the equation of motion for the particles in the ZA captures 
significant information that is encoded beyond linear theory. As we will argue, this can 
explain the large cross-correlation between the Zel'dovich density field and the true density 



field. 
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1.2 Beyond Standard Perturbation Theory 

Recently various analytical tricks have been put forward to remedy the convergence properties 
in the weakly nonlinear regime of SPT [8]. These include renormalized perturbation theory 
RPT [10], the path-integral approach [16], and the renormalization group flow [17]. Most 
of these expansion schemes work in the single-stream (or fluid) approximation. Thus, their 
results are applicable before shell-crossing, i.e. long before virialization can occur. They work 
only with the first two velocity moments of the Vlasov equation which reduce to the usual 
continuity and Euler equations. Since the higher moments of the one-particle distribution 
function are artificially discarded, one closes the system by introducing an equation of state, 
or equivalently a sound speed, which is set to zero for the CDM (along with any anisotropic 
stress) . The effects of the pressureless perfect fluid approximation on the power spectrum at 
the mildly non-linear scales relevant to the BAO at z = have been estimated to be around 
1% or more (e.g. [18-20]), which may not be sufficient to satisfy the requirements imposed 
by observations.^ 

Nevertheless, some of these expansion schemes have gained prominence and a discussion 
of their results is warranted. We will focus on two representative methods, one of them 
derived starting from the fluid approximation, while the other starting from LPT. Those are 
the renormalized perturbation theory (RPT) [10] and the Lagrangian Perturbation Theory 
version of [21].^ In these approaches (and most other recently proposed methods), one writes 
the non-linear power spectrum, Pnl-, as a sum of two pieces Pnl = R^Pl + Pmc-, where 
Pl is the power spectrum from linear theory, R is the response function (or "propagator") 
which gives the memory of the initial conditions, and Pmc gives the power generated from 
mode-coupling. This split has been extensively used recently, and has been applied in the 
reconstruction of the acoustic peak [4, 5], and so deserves special attention. 

In the split above, the response function is simply defined as the projection of the 
non-linear density field on the linear density field, and is therefore proportional to the cross- 
correlation between the initial density field (which is in turn proportional to the linear density 
field) and the density field at the moment in time of interest. However, as we already discussed 
that cross-correlation decays quickly in the mildly non- linear regime {pNL,L in Figure 2). 
Thus, R must decay quickly at small scales as well.^ This can be clearly seen in Figure 1 
where we plot R^Pl- 

Written that way, the mode-coupling power, Pmc-, has three sources: bulk flows, shell 
crossing and structure formation, where (for brevity) we define the latter to refer to all effects 
not captured by the ZA (see below). So, Pmc must compensate for the difference in power 
denoted by the green (long dashed) arrow in Figure 1. 

Let us try to understand the above split when written for the density field calculated in 
the ZA. The power spectrum in the ZA, P^, can be similarly split into Pz = R^Pl + Rmc 
where a Z indicates that we are working in the ZA; and Rz is approximately equal to R. 
The ZA leads to a similarly decaying density response function as for the non-linear power 
spectrum. However, the decay seen in the term R^Pl is completely absent in the total 
Zel'dovich power spectrum, Pz (see Figure 1), being compensated by the mode-coupling 

'^The effect of these corrections on the acoustic peak has not yet been evaluated, especially because all 
relevant calculations of the effect have been done using severe approximations. 

^To distinguish the original LPT from its approximated version offered by [21], we will refer to the latter 
as MLPT. 

®At high k, B? is well approximated by exp(— fc^cr^) [22], where (t„ ~ 9Mpc at 2: = is (1/^3) times the 
rms displacement of the particles in linear theory. 
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power P^^Q which is due to the effects of the large-scale flows and free streaming. This 
result can be understood in light of our discussion of the effects of the bulk flows: Although 
the memory of the initial conditions for the density at small scales is indeed lost under the 
random large-scale flows, the small-scale density correlations evaluated at the same moment 
in time are in fact preserved. Indeed, in the ZA one can show (Section 2) that the split 
above is due to the ability to choose different ordering parameters, over which to expand Pz- 
The expansion parameter resulting in the split of Pz above is in fact not small, since it is 
~ cr'^k'^ ~ 1 at the scales relevant for the BAO, where ~ 9Mpc at z = is related to the 
rms displacement of the particles in linear theory. So, the splitting of Pz into R^Pl + Pmc 
is not perturbative for the mildly non-linear regime, and this argument can analogously be 
extended to the splitting of P/vL in RPT (see also Section 2). To put it in another way, 
RPT and MLPT split the effect of the bulk flows, parametrized by kay{< k), between the 
pieces R^Pl and Pmc- However, as discussed above, the effects of the large-scale flows must 
cancel for equal-time statistics ~ something which is no longer true order by order in RPT 
and MLPT, although it must hold non-perturbatively. 

Thus, by using the density response function, both RPT and MLPT treat the cancella- 
tion of the effects of the bulk flows as something to be recovered only at higher orders. This 
problem can be circumvented by realizing that information about the large-scale flows (as 
well as some of the effect of shell crossing^^) can be easily restored by working in the ZA. 

In this paper we propose a different, physically motivated way to split the effects of 
mode-coupling, which will allow us to construct a simple model for the mildly non-linear 
regime. If one is able to collect the mode-coupling power due to the bulk flows with the 
power due to the memory of the initial conditions, R^Pl, one may be able to write Pnl = 
R^Pz + PmCj with a new R and new Pmc which are due to structure formation (by our 
definition). Here we defined R as the projection of the non- linear field on the Zel'dovich 
density field. Since the ZA captures the effects of the bulk flows^^, and R is proportional to 
Pnl,z (the factor of proportionality depending on equal-time statistics), we can conclude that 
all three quantities, R, Pz and Pmc, are independent of the large-scale fiows. Indeed, the 
Zel'dovich density field is well-correlated with the non-linear density field well into the mildly 
non- linear regime (see Figure 2), and therefore, the two are simple /c-dependent rescalings 
of each other. Thus, in Figure 1 one can clearly see that R?Pz follows Pnl closely even for 
k ^ k]\fL- The new Paic must be such as to compensate for the difference in power indicated 
by the blue (dot-dashed) arrow in Figure 1. So, we find that Pmc is nearly irrelevant for the 
BAO (Figure 1). Judging from the smaller mode-coupling power that needs to be obtained, 
calculating Pnl in the mildly non- linear regime using the split we proposed, may be easier 
than the route taken by RPT. 

In Section 2 we further discuss the proposed split of Pnl- In Section 3 we demonstrate 
how the proposed split can be applied in correcting N-body simulations for sample variance 
in the linear and mildly non- linear regimes. We show that the results presented in this 
paper will allow for a speed-up of an order of magnitude (or more) in the scanning of the 
cosmological parameter space for studies which require a good handle of the mildly non- 

^°The equation of motion corresponding to LPT treats the density in the Poisson equation in the single- 
stream approximation. However, the resulting particle trajectories are still allowed to intersect, resulting in 
a rich structure in phase-space. It is in this sense that we say that LPT, and the ZA in particular, capture 
some of the shell-crossing information but not all of it. 

^^As shown by [23], the ZA satisfies Galilean invariance for the continuity equation and the divergence of 
the Euler equation. 
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linear regime, such as those targeting the BAO. In Section 4 we explain why the method 
we propose for eliminating the sample variance from simulations works so well. Then in 
Section 5 we construct a simple physically motivated model for the mildly- nonlinear regime, 
which captures the effects of bulk flows^^; and which models the matter power spectrum 
to 1% accuracy for these scales and gives robust estimates for the acoustic scale. We then 
summarize our results in Section 6. In Appendix A we provide some more details on the 
derivation of our model. 



2 Splitting the non-linear power spectrum 

As we discussed in Section 1 , many recently proposed expansion schemes targeting the mildly 
non-linear regime, such as the RPT and MLPT expand the non-linear power spectrum, Pnl, 
as a sum of two terms - a term proportional to the linear power spectrum, Pl, and a 
remainder, called the mode-coupling term, Pmc- 

PNL{k) = RHk,rj)PL{k,7])+PMc{k,ll) , (2.1) 

where rj is conformal time. The response function, R is given by the projection of the non- 
linear density field, Jatl, on the linear density field, Sl- 

= (2.2) 

Clearly, R can be obtained directly from the cross-correlation coefficient, pnl,l which is 
plotted in Figure 2: 

^^(^,^) = ^^pkL(^,^)- (2.3) 

To understand the above split better, let us write it for the density field in the Zel'dovich 
approximation (ZA) as we did in Section 1. The Zel'dovich power spectrum, Pz can be 
written as: 

Pz = RlPL + Piw , (2.4) 

with Rz defined in analogy with R in eq. (2.2) but with Snl replaced by the Zel'dovich 
density field, 6z- The exact expression for the matter power spectrum, Pz{k), in the ZA is: 

{6z{k)6z{k')) = 6D{k + k')Pz{k) (2.5) 
= 6D{k + k') I ^e-*''exp[-A:,A;,(5,,a2-vI/,^.(g))] , 

where 

I d'^^-^PLiOe''-^ (2.6) 

is proportional to the linear velocity correlation function. The velocity dispersion is propor- 
tional to (T„, which is defined as = ^'m(0)/3. 



In the sense that bulk flows are irrelevant for equal time-statistics order by order in the proposed split of 
Pnl- 
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At first order in P^^ we liave Pz ~ Pl as expected. At second order, one can check 
explicitly that the large-scale modes affect Pz only through the density power spectrum, and 
not through the velocity power spectrum. Thus, the ZA is not affected by bulk flows at 
second order. This result can be extended to all orders in the overdensity [23]. 

If we follow [21] to recover the MLTP result, we should expand Pz in ^ij, while keeping 
(j^fc^ in the exponent. The first order result in this expansion becomes exp(— /c^o"^)^^ which 
we can identify with R^Pl [24]. The rest of the terms must be collected in Pmc- 

Therefore, the split above distributes the canceling contributions from the bulk flows 
(terms containing cr^(< k)) between R\Pl and P^c (see also Section 1). The expansion in 
^ above is in the parameter ~ fc^o"^ which is ~ 1 in the mildly non-linear regime. Thus, we 
can expect that P^q comparable to if not larger than R^Pl, which is indeed the case (see 
Figure 1, keeping in mind that Rz and R are almost equal). 

One can apply the analysis above to the non-linear density field as well, noting that Rz 
is a very close approximation to R [22]. Thus, we can conclude that the split, eq. (2.1), is not 
optimal for analyzing -P/vL- Therefore, for example, by choosing an initial power spectrum 
with spectral index —3 < n < —1, one can arrange for ct^(< k)k'^ to diverge, while the 
overdensity remains finite. In that limit, RPT and MLPT break down completely; while the 
split proposed below avoids these divergences, as equal-time statistics cannot depend on the 
large-scale bulk flows. 

So, instead of writing the above split around Pl, we perturb around the Zel'dovich 
approximation, which captures the effects of the bulk flows (see Section 1; and [23] as well). 
Thus, we will keep Pz exact, and expand Pnl around it: 

Pnl = R\k, 7j)Pz{k, rj) + PA4c{k, v) , (2-7) 

with 

Rik,,)^%§^. (2.8) 

In [25] it was shown that the splitting above emerges naturally when considering the statistical 
properties of CDM in phase-space. 



3 Correcting for sample variance in N-body simulations 

The Zel'dovich density field and the non-linear density field are very well cross-correlated 
across the mildly-nonlinear regime (Figure 2). Thus, one can treat them as simple rescalings 
of each other, the rescaling given by R. 

As we discussed in Section 1, we expect Pz to capture most of the effects of the large- 
scale modes, which are the ones most affected by sample variance. Thus, let us assume 
for now that R is not affected by sample variance, and therefore does not depend on the 
particular realization of the random density field (see Section 4). We can then write: 

SNL{k, 7]) R{k, r])6z{k, rj) + 6Mc{k, v) , (3-1) 

where 5nl and 6z are calculated in the same realization, and Smc is defined such that 

{SmcS*z) = . 
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Thus, for each simulation box of S^iik,?]), one has to run a corresponding Zel'dovich "sim- 
ulation" with the same initial conditions to get the corresponding 6z{k,ri). 

Note that eq. (3.1) can be treated as a linear model for 5]srL{Sz), with a non-random 
coefficient R{k,r]) and a random residual given by Smc- The least-squares fit gives the 
following estimators (denoted with hats) for R and Pmc- 



R 



MC 



R\Sz5*z) 



Thus, we propose the following estimator for Pnl- 



pNL{k,rj)=R'Pz + P, 



MC 



{^nlS*z] 
{SzS*z) 



Pz - {SzS*z) + {SnlS[ 



(3.2) 
(3.3) 

(3.4) 



where Pz is defined below (eq. (3.5)). The averages above denoted by the angular brackets 
are averages both over the modes in each A;-shell in each realization, and over all simulation 
boxes. Therefore, the estimators above are manifestly unbiased in the limit of an infinite 
number of (finite in size) simulation boxes. 

The Zel'dovich power spectrum Pz on the right-hand side of eq. (3.4) is the exact 
theoretical power spectrum, Pj (we put a subscript T to distinguish the power spectrum 
given in eq. (2.5) from the Nyquist-corrected Pz below), corrected for the truncation of Pl 
at the Nyquist wave-vector, k^. A good estimate for Pz is given by: 



Pzik,v) ~ Pzik,v) exp 



47r , 



V 



Pl{w, ri)dw 



(3.5) 



which corresponds to (2.5) with a velocity dispersion, o"^, evaluated up to k^- Alternatively, 
one can compute Pz from an ensemble of initial conditions, since each Zel'dovich "simulation" 
is cheap to compute. The latter method is slower, but captures the effects of the cubical cells 
for initial conditions on a grid, or the effects from glass-like initial conditions. 

By applying eq. (3.4) we will show that we can extract P^Lik, rf) with a great accuracy 
even from a single N-body simulation, as long as we calculate the Zel'dovich density field 
using the same initial conditions used for the fully non-linear N-body simulation. As we will 
show below, one can even do better by applying eq. (3.4) using the overdensity field obtained 
at second order in LPT (2LPT), instead of the Zel'dovich overdensity: 



PNL{k,vi) 



{^nlSIlpt) 
{62lptS2lpt) 



P2LPT — {^2LPt52LPt) + {^NlS 



-'NL/ 



(3.6) 



The exact P2LPT above can be computed from a large ensemble of 2LPT "simulations", 
which are cheap to compute. 

In order to test the estimators for Pnl we ran 10 simulations using GADGET-2 [26] 
with 256^ particles each, in a box with side Lbox = 500 /i/Mpc. The cosmological parameters 
we used are (in the standard notation) as follows 



(1^6,J7matter,J^A,/i,n„cJ8) = (0.046,0.28,0.72,0.70,0.96,0.82) . 



(3.7) 
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Figure 4: We show the raw non-hnear matter power spectrum at z = (with 2-sigma errorbars for 
the average Pnl) obtained from an ensemble of 10 N-body simulations. We also show the linear and 
Zel'dovich power spectra for comparison. All power spectra are divided by a smooth BBKS power 
spectrum to highlight the BAO wiggles. 

The initial conditions were set at a redshift of 49, using the 2LPT code provided by [27]. 
The density field was calculated on a 512^ grid using Cloud-in-Cell (CIC) assignments, and 
the final power spectra were CIC corrected. 

In Figure 4 we show the non-linear power spectrum calculated directly from the simu- 
lations. Note the large errorbars on the average of Pnl which are a direct consequence of 
sample variance. In Figure 5 we plot the same quantities as in Figure 4 obtained from the 
same 10 simulations but corrected using the 2LPT scheme, eq. (3.6). To do that, we used the 
following method: For each of the 10 N-body simulations, we ran a 2LPT "simulation" with 
the same initial conditions. To obtain Pnl we used the 2LPT estimator, eq. (3.6), with the 
averages running over both all /c-modes and all simulation boxes. The true P2LPT required 
for the 2LPT estimator was obtained from an average of 400 2LPT simulations. Then to cal- 
culate the errorbars, we calculated Pnl using eq. (3.6) from each of the 10 realizations; took 
the rms of the resulting ten Pnl estimates, and divided by -v/lO to obtain the 1-sigma error 
in the mean Pnl- The results are shown in Figure 5. One can clearly see that the sample 
variance has been practically eliminated in the linear regime; and is significantly reduced in 
the mildly non-linear regime. 

In Figure 6 we show the amount of improvement in the errors in determining Pn l using 
eq. (3.4) applied with 2LPT (i.e. eq. (3.6)), the Zel'dovich approximation, and linear theory 
(by replacing the quantities evaluated in the ZA in eq. (3.4) by their counterparts in 2LPT 
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Figure 5: The sanie as in Figure 4 with the non-linear power spectrum calculated from the same 
10 simulations but using the 2LPT estimator, eq. (3.6). 



and linear theory). In the BAO region (0.05/i/Mpc< k < 0.15/i/Mpc) the 2LPT estimator 
results in a reduction of the errorbars by about 5-50 times (depending on the scale) compared 
to the standard result. This implies that to achieve the same errors in the power spectrum in 
that region, the necessary computational time is reduced by a factor of ~ 30 — 2000. Linear 
theory produces much worse improvement than LPT due to the fact that the cross-correlation 
Pnl,l deviates from 1 at much larger scales (Figure 2; see also Section 4). As discussed in 
Section 1, this in turn can be understood by the fact that unlike LPT, linear theory does not 
properly account for the random bulk motions, which decorrelate the linear and non-linear 
density fields. 

Note that the improvements obtained in both 2LPT and the ZA are no longer significant 
above k ~ 0.5/i/Mpc. This is the scale where both cross-correlations pnl,z and pnl,2LPT 
become ~ 0.9 (see Figure 2 for pnl,z)- Neither the ZA nor 2LPT can capture the small- 
scale dynamics correctly. Thus, Pmc grows large beyond k ~ 0.5/i/Mpc, and the Pnl 
estimators we proposed no longer result in error improvements. However, those scales are 
usually well within the simulation volumes, and thus sample variance is usually less of a 
problem. Moreover, those scales are irrelevant for the acoustic peak from the BAO, since the 
oscillations in the power spectrum are completely washed away at these scales. 
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Figure 6: We show the improvement in the error, a, in determining Pnl &t z — using the restora- 
tion schemes discussed in the text: using eq. (3.4) apphed with 2LPT, the Zel'dovich approximation, 
and hnear theory. We denote the error in Pnl as obtained directly from the N-body simulations as 

4 Why does the method work? 

We need to understand why the procedure described in the previous section works so well in 
accounting for sample variance. For example, one may wonder why linear theory works so 
much worse than the ZA, when it models Pnl even better than Pz (Figure 1). Moreover, 
we have to try to understand why 2LPT works so much better than the ZA. 

As long as the cross-correlation coefficient ppt,nl is close to 1 (here PT stands for 2LPT, 
Zel'dovich or linear theory), 6nl and 5pt are simple rescalings of each other. Therefore, we 
can write Snl ^ R6pt (where for brevity we drop the tilde from R and Pmc whenever used in 
conjunction with the subscript PT, irrespective of whether PT stands for LPT or not) with a 
non-random R. Our P^l estimator boils down to Pnl ~ R^PpT with R^ {\^nl\'^) / {\Spt\'^) ■ 
Thus, the Pnl estimator has vanishing sample variance since for ppt,nl ^ 1, R must be 
non-random; and PpT is non-random by construction. 

The above argument implies that when ppt,nl ~ 1 we must have a negligible Pmc 0,^(1 
the fractional error due to sample variance of R?Ppt (which equals that of R?, since Ppt is 
the theoretical PT power spectrum) must be sub-Gaussian, i.e. less than y^2/Nk, where 
is the number of modes in each fe-bin from all simulations. To test the first prediction, in 
Figure 7 we show the ratio of the two terms entering in the split of Pnl, eq. (2.7), for the 
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Figure 7: We show the ratio of the two contributions to Pnli eq- (2.7), depending on the choice of 
split (around the linear, Zel'dovich or 2LPT power spectrum) of Pnl- Note that the mode-coupling 
power is negligible at large scales, and becomes important only for scales for which ppt,nl deviates 
significantly from one. 



split applied around Pl, Pz and P2LPT- One can see that at large scales, while ppt,nl ~ Ij 
it is indeed B?Ppt that dominates over Pmc- To test the second prediction, in Figure 8 we 
plot the fractional errors of and Pmc- As expected, we see that as long as ppt,nl ~ 1 
the fractional error of B?' is much smaller than the variance for a Gaussian random variable, 
\/2/Nk. This is not the case for Pmc which has slightly super-Gaussian errors. However, 
Pmc is negligible at large scales. 

We can conclude that for scales with ppt,nl ~ 1 the overall error of Pnl must be 
greatly suppressed, which is indeed what we saw in the previous section. Moreover, since it 
is the large-scale modes which are most affected by sample variance, one can conclude that 
it is R and PpT, and not PmC: that carry the information on the coupling to the large-scale 
modes. We will elaborate on this point further below. 

As we saw, having ppt,nl ~ 1 is crucial for the proposed split of Pnl-, eq. (2.7), to 
work well. We can restate that condition as follows. An adequate reconstruction method 
for Pnl should rely on a split of Pnl such that: 1) the mode-coupling power is small in the 
mildly non- linear regime, which implies that the approximate density field around which we 
expand (in the previous section that was the field in the ZA or 2LPT) is well correlated with 
the non-linear density field; 2) the power spectrum PpT of the approximate density field does 
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Figure 8: We show the fractional error of i?^ and Pmc (symbols), as well as the naive estimate of 
the error for Gaussian random variables, given by ^J2/Nk (solid line). Note that the error of R is 
greatly suppressed for ppt,nl ~ 1- 

not deviate much from Pnl at the scales which have large sample variance. Both conditions 
above are satisfied in the ZA and 2LPT but not in linear theory. Linear theory has a large 
Pmc in the mildly non-linear regime, which is a direct consequence of the pl,nl not being 
close to one. This explains why the split around Pl is inferior to a split around the Zel'dovich 
or 2LPT power spectra. 

Next, let us try to understand why 2LPT works so much better than the ZA. To do 
that we performed a series of numerical experiments. First we chose a scale kc = 0.15 /i/Mpc 
which is accessible using perturbation theory. We then ran a set of 10 simulations in which 
we cut off the high-k modes from the initial conditions (IC) by setting Pinitiai(^ > fee) = 
- we will refer to this set as the "long" set; and another set of 10 simulations in which we 
removed all large-scale modes from the IC: -Pinitiai(^ < fee) = - we will refer to this set 
as the "short" set. We correct both sets using eq. (3.6) with the S2LPT calculated with the 
correspondingly band-passed initial conditions. 

Mode coupling generates power even at those scales, which had zero initial power. In 
Figure 9 we show the fraction of the full Pnl which comes from the long and short IC, as well 
as the residual cross power, which arises from the coupling of the initially present short-scale 
and large-scale modes: 

PNL,iun = Pn L,long + Pn L,shovt + -P/VL,cross , (4.1) 

where a subscript "full" implies the quantity evaluated with the full (both large- and short- 
scale modes included) initial conditions, i.e. PArL,fuii = Pnl- 

As one can see from Figure 9 practically no large-scale power is generated in the 
simulations with initially only short-scale modes. The reason is that mass and momen- 
tum conservation require that any power leaking from the small scales to the large scales 
should scale at least as (/co"^(> kc))'^ [28], where av{> kc) ~ 2.8Mpc//i is the rms par- 
ticle displacement due to the small scale {k > kc) power. This scaling is modified to 
~ {k/kc)'^{kav{> kc))'^ ~ {k/kc)^/32 in Figure 9 since we have divided by Pnl which scales 
as roughly k~'^ around kc- This results in the sharp cutoff observed for PNL,short for k < kc- 
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Figure 9: We show the fraction of the fuU Pnl at z = generated from the long and short IC, and 
the residual cross power spectrum. Errorbars are suppressed for clarity, but all curves are corrected 
for sample variance using eq. (3.6). The position of the cutoff kc is denoted by a vertical line. 



The cross power generated from the coupling of the initially long and short modes is 
again mostly localized at high k > kc- This term tells us how the small-scales are influenced 
by the effects of the perturbed local J^matter and by the tidal (and higher order) interactions 
due to the presence of the large-scale modes. 

The simulation with low-pass filtered IC produces a fat tail at high k > kc- This is a 
direct consequence of mode-coupling both due to structure formation and stream crossing. 

For k < kc we can see that it is the low-pass filtered IC that dominate the true power 
spectrum, Pnl- Thus, we can conclude that the power in the mildly non- linear regime not 
captured by linear theory is due almost entirely to the coupling to the large-scale modes. 

Since the ZA and 2LPT satisfy both requirements above for an adequate reconstruction 
method, they should account for much of the large-scale power in Pnl, as well as the peak 
in the cross power at kc, which at these scales should mostly be due to large-scale flows. To 
investigate this further, in Figures 10 and 11 we plot the different contributions, i, to the 
full non- linear power which are not accounted for by the Zel'dovich/2LPT approximation: 
{PNL,i - Pz/2LPT,i)/PNL,fun- ^or negligible Pmc this quantity reduces to 

PNL,i — Pz/2LPT,i ^ Rf — 1 Pz/2LPT,i 
PNL,iu\\ -Rf^jj Pz/2LPTM^ 
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Figure 10: We show the fractional contributions to the fuU P^l at z = which are not accounted 
for by the ZA: those contributions generated from the long and short IC, and the residual cross power 
spectrum. Errorbars are suppressed for clarity, but all curves are corrected for sample variance using 
eq. (3.6). The position of the cutoff is denoted by a vertical line. 

where a subscript "full" implies the quantity evaluated with the full initial conditions as 
before; and R is defined either with respect to the ZA or 2LPT. 

Comparing Figures 10 and 11 to Figure 9 we can see that the peak of the cross power 
at /c ~ A^c is well-accounted for by LPT (both the ZA and 2LPT). So, indeed LPT captures 
quite well the cross-correlations between nearby large-scale modes. The residual — 1 is due 
mostly to large-scale modes in the ZA (for k < kc), while P2LPT captures better the coupling 
between the large-scale modes. Since it is precisely the large-scale modes which suffer most 
from sample variance, this explains why 2LPT is much better than the ZA in correcting for 
sample variance (see Figure 6). 

5 Modeling the transfer function in the mildly non-linear regime 

In the previous sections we demonstrated the advantages of splitting the non-linear power 
spectrum as in eq. (2.7). Since the mode-coupling term, Pmc^ is relatively small in the 
mildly non-linear regime, we can think of i?(/c, rj) as a transfer function transforming the 
power spectrum in the ZA to the true Pml- In this section we present a simple physically 
motivated model for the mildly non-linear regime from which we obtain a good approximation 
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Figure 11: As in Figure 10 but for the the fractional contributions to the fuU Pnl which are not 
accounted for by 2LPT. 

for R. In this section we will work only with R as written in eq. (2.7), i.e. with respect to 
the ZA, not 2LPT. 

From Figure 10 and eq. (4.2) one can see that ii^ — 1 is dominated by large-scale 
modes. Those affect the smaller scales by rescaling the local matter density (i.e. producing 
perturbations to the local r^matter); through tidal effects; and through higher order (in a 
multipole expansion) corrections. We will concentrate on the leading corrections, which 
correspond to those coming from the perturbed local r^mattcr- These perturbations lead to 
a perturbed local growth factor, whose effect is analyzed in [29]. These authors split the 
density field {5 = 5\ + 5s) into a large-scale ((5a) and a short-scale {5s) component and find 
that the large-scale modes modify the short-scale modes as follows (in configuration space): 

5s,NL « 5s,L [l + -5k,l + y^'^IlJ (5-1) 

plus corrections which are higher order in 5j^^l (for /cshort ^ ^long, where the A;'s correspond 
to 5s and 5\, respectively). As usual above we used a subscript L to denote quantities in 
linear theory, while NL denotes the "true" result. The rescaling above due to the large-scale 
modes comes exactly from the corrections to the local growth factor from the slightly changed 
^matter- Assuming a sharp fc-cutoff for splitting 5 into 5s and 5\, the resulting short-scale 
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Figure 12: Here we show various matter power spectra at z = for ACDM. As before, the power 
spectra are divided by a smooth BBKS power spectrum. The non-hnear power spectrum is given 
by P/vl; the hnear power spectrum by P^; the 1-loop result in SPT is denoted by ISPT; the power 
spectrum calculated in the Zel'dovich approximation is Pz- Our formula for the power spectrum 
in the mildly non- linear regime results in a power spectrum given by the curve denoted with B?Pz 
(eq. (5.5)), with ci = 2950/1323 ~ 2.23. The shaded region is bounded with the Pnl^ obtained from 
(5.3) with A = 0.9fc/2 (lower bound) and A = l.lfc/2 (upper bound). 



power spectrum is 

Pnl « + 4^ / PL{k)k'dk\i = (^1 + II 5\< A)) , (5.2) 

where we used that A:short ^ ^long! A is a large-scale cutoff; and the linear short-scale and 
large-scale modes are independent: {6s,l^a,l) = 0. As before (5^(< A) is the density variance 
integrated up to A. 

In order to recover R, we should rewrite the above equation using Pz, where for brevity 
we drop the subscript T from Pj. However, the ZA already captures some of the effects 
from (5a- By expanding Pz entering in Pnl ~ R^Pz in powers of 6'^{< A) and matching the 
resulting Pnl to the short-scale power in eq. (5.2) we obtain (see Appendix A) 

PNLik, 11) « Pz{k, ^) (l + II '^'(< A)) . (5.3) 

Note that the above equation is derived from the following expression relating the short- 
scale non-linear overdensity and the overdensity in the ZA valid for equal-time statistics (see 
Appendix A): 

8s,NL = 5z{l + hK,L + ^5lA ■ (5.4) 
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This equation is our model for the overdensity in the mildly non-linear regime. It is correct 
to second order in 8\^l (for ^short ^ Mong)- However, when using the above expression one 
needs to keep in mind that {5z6a,l) 7^ 0, unlike {6s,lSa,l) = 0. 

The above results surely break down for A > k/2, since then two large-scale modes can 
couple to directly produce a short-scale mode. Thus, we choose A = k/2: 



PNL{k,r,) - Pz{k,r,) 



(5.5) 



from which one can read off E? as the term in the brackets. We plot this result in Figure 12. 
For reference we also give the corresponding 2-pt functions in real space in Figure 13. 
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Figure 13: Here we show the 2-pt function, ^, at z = for ACDM. One can see the acoustic peak at 
« 105Mpc//i. The dotted curve is the ^ corresponding to R^Pz with R given by our model, eq. (5.5). 

From Figure 12 one can see that the expression for the matter power spectrum, eq. (5.5), 
reproduces the true power spectrum to 1% accuracy for k < 0.2/i/Mpc. Thus our simple 
model is working much better than 1-loop SPT and has about the same level of accuracy as 
2-loop RPT in the mildly non-linear regime [6]. 

One may argue that the (1/2) in (5^(< k/2) in eq. (5.5) should really be treated as a 
fudge factor. Since in real space the acoustic peak is shifted by about 1% from the result in 
linear theory, we varied the coefficient of (1/2) to produce a comparable shift. We found that 
deviations by as much as 50% are needed. ^"^ These result in wildly different power spectra (see 
Fig. 12 where with the shaded region we show the result from a 10% change in the coefficient 
of (1/2)). Thus, barring any systematics, the model behind eq. (5.5) gives a robust estimate 
for the acoustic scale and the mildly non- linear power spectrum. Furthermore, one may be 
able to apply our model for the reconstruction of the BAO peak in an analogous way to [4] . 

Before we conclude, a word of caution is in order. In order to simplify the calculations 
in this paper, we assumed that the large-scale overdensity is spherically symmetric, eq. (A. 4), 

^^Varying the coefficient in front of <5^(< k/2) by 50% results in even smaller shifts in the acoustic scale. 
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and not a plane-wave as it should be. This assumption enters in the derivation of both 5z, 
eq. (A. 6), and Ss,nl, eq. (5.1). In principle, this assumption can lead to problems. And indeed 
eq. (A. 7) tells us that Pz > Pl, which is clearly not the case (e.g. Fig. 13). However, we 
use Pl only for convenience in the intermediate results, hoping that any mistake we make, 
affects Pz and Pnl equally (because of the strong correlation between 6z and Snl)i and 
therefore may cancel in the final result, (5.3,5.5). Indeed our numerical results confirm this 
intuition. Clearly, if one wants to do better, one can repeat the calculation with a plane-wave 
large-scale mode. However, we leave that for future work, noting that our analytical formula 
(5.5) gives rather robust results as it is. 

Therefore, we can conclude that the CDM particle dynamics in the mildly non-linear 
regime is well approximated by the dynamics in the Zel'dovich approximation with perturbed 
growth factor. These perturbations to the local growth factor result from the perturbed local 
matter density caused by the presence of the large-scale modes. 

6 Summary 

In this paper we used LPT to investigate the behavior of modes in the mildly non-linear 
regime. We find that their behavior is affected mainly by neighboring (in k) modes and 
by larger-scale modes. Much-smaller-scale modes are isolated and have almost negligible 
back-reaction. At the same time large-scale modes affect the mildly non-linear scales only by 
changing the local matter density or through subdominant tidal effects. 

We find that the effects of the neighboring modes and large-scale modes are well modeled 
by LPT (but see below), which takes into account the effects of the bulk flows. Thus, 
unlike linear theory, the cross-correlation between the LPT overdensity and the non-linear 
overdensity is close to one in the mildly non- linear regime. We also find that the ratio (which 
approximately gives the square of the response (transfer) function, R) of the true power 
spectrum and the power spectrum in the ZA is nearly independent of realizations, which is 
due to the fact that the behavior of the large-scale modes, which are the ones with the large 
sample variance, is well-captured by LPT in the mildly non-linear regime. 

Using the above results, we developed a straightforward method for correcting the results 
from N-body simulations for sample variance in the linear and mildly non-linear regimes. The 
method has no free parameters and introduces about an order of magnitude improvement in 
the errors in determining the true matter power spectrum. The method can be used to speed 
up the scanning of the cosmological parameter space by an order of magnitude or more for 
observables in the (linear and) mildly non-linear regime. 

We then constructed a physically motivated model for the mildly non-linear regime, 
which treats the effects of the bulk fiows correctly. The model is the same as the Zel'dovich 
approximation but with a perturbed local growth factor. Those perturbations are caused 
by the large-scale modes changing the local ^matter - an effect which (to first order in the 
large-scale power) is not completely captured by the ZA and 2LPT. 

Using this model, we gave an approximate expression for the transfer function, R, which 
relates the power spectrum in the ZA, and the true power spectrum. The resulting power 
spectrum models the true power spectrum with an accuracy of < 1% in the mildly non-linear 
regime, and gives robust estimates for the BAO scale. 
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A Deriving equation (5.3) 



In this section we derive equation (5.3). To do that let us start by writing down the expression 
for the Eulerian comoving coordinates x{q,7]) of a particle in the ZA, which are a function 
of the particle's Lagrangian coordinates q and rj: 

x{q,v) = q + D{r])s{q) , (A.l) 

where D is the growth factor, and s is an irrotational displacement field given by 

-D{r,)Vg-s{q)=6L{q,v) ■ (A.2) 
The overdensity in the ZA, is then given by {k ^ 0): 



5zik, r^) = j (^e— j d'qSo {x-q- D8{q)) = j ^e-'=-^^^^'^^+''^ . (A.3) 

We can write the displacement field as a sum of a short-scale, s^, and a large-scale piece, 
s\. To lowest order in A;iong, the large-scale displacement can be written as (in Lagrangian 
space): 

SK{q,r]) ~ ~^;3]^'^a,l(?/) , (A.4) 

which corresponds to a (5a, L being constant in space. Thus, in our calculation we only take 
into account the change in the local r^matter induced by the presence of (5a,l- 

Then we can expand 5z to third order in s (i.e. third order in 5l)^ and use the split 
s = Sg + Sf^ together with eq. (A.4). After some algebra we obtain (in configuration space) 

5z = 5,,l{i + \6k,l+^^51j\ (A.5) 



5 L>2 , 

+D (sa • Vq) 5s,L - -^^^KL (SA • Vq) 5s,L + ^S^SA^g^ 9gi (5s,L 



Now let us show that the second line gives no contribution to the power spectrum Pz cal- 
culated to second order in P^. When calculating Pz, we use that the linear short-scale and 
large-scale modes are independent ((<5s,l(5a,l)=0). Thus, we can see that the second term on 
the second line gives a contribution to Pz proportional to 

{5K,Ls\){6s,Ldqi6s,L) = , 

which vanishes because of the vanishing first bracket (by isotropy). 

Analogously, the first and third term on the second line of eq. (A.5) give equal but 
opposite in sign contributions to the power spectrum, which therefore cancel out as well. 
Thus, we can see that at this order Pz depends only on 5a,l but not on the large-scale 
displacement field. This must indeed be the case since as we already argued equal-time 
statistics in the ZA must be independent of the large-scale bulk fiows. 

Thus, for the purposes of equal time statistics, we can simply write 

= <5.,L (l + ^5a,l + y , (A.6) 
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plus irrelevant bulk flow terms. The power spectrum is then given by 

Pz{k,7]) = PL{k,rj){l + A6\<A,7])) . (A.7) 

Expressing Pl from the above equation and plugging it in eq. (5.2) we recover eq. (5.3), 
which was the goal of this section. Note that by comparing eq. (A. 6) and eq. (5.1) we obtain 
eq. (5.4), which is our model for the short-scale non-linear overdensity. 
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